      PROGRAM DRIVER
      USE AERO_CONFIG
      USE AERO_COAG, ONLY: TEST_COAG_COEF
      USE AERO_INIT  
      USE AERO_SETUP  
      USE AERO_SUBS
      USE AERO_NPF   
      USE AERO_DIAM
      USE AERO_DISCRETE, ONLY: DISCRETE, DISCRETE_OUT
      USE AERO_ACTV, ONLY: GETACTFRAC
      USE AMP_AEROSOL
      IMPLICIT NONE
      INTEGER :: I,J,N,IT,NIT,IREP,NREP,NSETS,ICSET,SET,ITMATRIX
      REAL(8) :: AERO(NAEROBOX)
      REAL(8) :: GAS(NGASES)
      REAL(8) :: NACT(NMODES)
      REAL(8) :: TK                    ! absolute temperature [K]
      REAL(8) :: RH                    ! relative humidity    [0-1]
      REAL(8) :: PRES                  ! ambient pressure     [Pa]  
      REAL(8) :: TSTEP                 ! model time step [s]
      REAL(8) :: TSUBSTEP              ! increment of model time step [s]
      REAL(8) :: EMIS_MASS(NEMIS_SPCS) ! mass emission rates [ug/m^3/s]
      REAL(8) :: AQSO4RATE             ! incloud SO4 production rate [ug/m^3/s]
      REAL(8) :: TSTOP                 ! stop  time [h]
      REAL(8) :: PAQSO4                ! incloud sulfate production rate during the time step [ugSO4/m^3/s]
      REAL(8) :: EMIS_MASS_IN          ! mass emissions rate for each time step (all species) [ug/m^3/s]
      REAL(8) :: PH2SO4,PHNO3,PNH3     ! production rates [ug/m^3/s]
      REAL(8) :: MH2SO4,MHNO3,MNH3     ! initial mass concentrations [ug/m^3]
      REAL(8) :: SPCMASS(NMASS_SPCS+2) ! total mass conc. of each model species [ug/m^3]
      REAL(8) :: T                     ! model time [h]
      REAL(8) :: DT                    ! model time step [h] ( = TSTEP / 3600.0D+00 )
      REAL(8) :: WUPDRAFT = 0.5D+00    ! cloud updraft velocity [m/s]
      REAL(8) :: DIAG(NDIAG_AERO,NAEROBOX)        ! budget or tendency diagnostics [ug/m^3/s] or [#/m^3/s]
      REAL(8) :: DIAG_ACCUM(NDIAG_AERO,NAEROBOX)  ! budget or tendency diagnostics [ug/m^3/s] or [#/m^3/s]
      REAL(8) :: SUMTMP                           ! scratch variable for DIAG_ACCUM sum
      REAL(8), PARAMETER :: TEMIS_MAX = 24.0d+00  ! turn-off time for emissions [h]
      LOGICAL, PARAMETER :: WRITE_OUT = .TRUE.
      LOGICAL, PARAMETER :: DO_TEST_COAG  = .FALSE.
      LOGICAL, PARAMETER :: DO_TEST_QUAD  = .FALSE.
      LOGICAL, PARAMETER :: DO_TEST_ACTV  = .FALSE.
      LOGICAL, PARAMETER :: DO_TEST_DEPV  = .FALSE.
      LOGICAL, PARAMETER :: DO_EVAL_DEPV  = .FALSE.
      LOGICAL, PARAMETER :: DO_DIAG_WRITE = .TRUE. 
      CHARACTER(LEN= 1) :: DUMCHAR
      CHARACTER(LEN= 1) :: QUOTE = "'"
      CHARACTER(LEN=80) :: FILENAME
      CHARACTER(LEN= 2) :: NSETCHAR
      CHARACTER(LEN= 2) :: DATCONV
      CHARACTER(LEN= 5) :: MECHCHAR
      
      INTEGER, PARAMETER :: UNIT1  = 10
      INTEGER, PARAMETER :: UNIT2  =  9
      INTEGER, PARAMETER :: NFILES = 16
      CHARACTER(LEN=4), DIMENSION(NFILES) :: EXTCHAR = (/'mssa','mssb','nmba','nmbb',
     &                                                   'so4a','so4b','bcar','ocar',
     &                                                   'dust','seas','mdpa','mdpb',
     &                                                   'spcm','gass','ccna','ccnb'/)
      INTEGER, DIMENSION(NFILES) :: OUTUNIT = (/11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26/)

      INTERFACE
        SUBROUTINE OUT(AERO,GAS,NACT,SPCMASS,T,OUTUNIT)
          REAL(8) :: AERO(:)          
          REAL(8) :: GAS(:)           
          REAL(8) :: NACT(:)           
          REAL(8) :: SPCMASS(:)               
          REAL(8) :: T
          INTEGER :: OUTUNIT(:)
        END SUBROUTINE OUT
      END INTERFACE

      EMIS_MASS(:) = 0.0D+00
      IXXX = 1
      IYYY = 1
      IF( MECH .EQ. 1 ) MECHCHAR = 'mech1'
      IF( MECH .EQ. 2 ) MECHCHAR = 'mech2'
      IF( MECH .EQ. 3 ) MECHCHAR = 'mech3'
      IF( MECH .EQ. 4 ) MECHCHAR = 'mech4'
      IF( MECH .EQ. 5 ) MECHCHAR = 'mech5'
      IF( MECH .EQ. 6 ) MECHCHAR = 'mech6'
      IF( MECH .EQ. 7 ) MECHCHAR = 'mech7'
      IF( MECH .EQ. 8 ) MECHCHAR = 'mech8'

      IF ( CPU_STATS ) CALL CPUTIME( 0 )

      OPEN(UNIT1, STATUS='OLD',    FILE='aero_driver.in')
      OPEN(AUNIT1,STATUS='REPLACE',FILE='logmatrx.out')
      OPEN(31,    STATUS='REPLACE',FILE='massadj1.out')
      OPEN(32,    STATUS='REPLACE',FILE='massadj2.out')
      OPEN(33,    STATUS='REPLACE',FILE='massadj3.out')
      OPEN(34,    STATUS='REPLACE',FILE='npfmtest.out')
      FILENAME = MECHCHAR // 'masstest.out'
      OPEN(35,    STATUS='REPLACE',FILE=TRIM(FILENAME))
      FILENAME = MECHCHAR // 'actvtest.out'
      OPEN(40,    STATUS='REPLACE',FILE=TRIM(FILENAME))
      OPEN(41,    STATUS='REPLACE',FILE='depvtest.out')
      OPEN(44,    STATUS='REPLACE',FILE='logdiags.out')
      IF( DISCRETE_EVAL_OPTION ) OPEN(36,    STATUS='REPLACE',FILE='logdiscr.out')
!     OPEN(UNIT2, STATUS='REPLACE',FILE='diffscript')

      READ(UNIT1,*) NSETS   
      READ(UNIT1,*) DUMCHAR

      CALL SETUP_CONFIG       
      CALL SETUP_SPECIES_MAPS 
      CALL SETUP_DP0           ! need default particle diameters for plots
      CALL SETUP_DIAM          ! need default particle diameters for plots

      DO N=1, NSETS

        READ ( UNIT1,*) SET,ICSET,ILAY,NREP,TSTEP,TSTOP,TK,RH,PRES,EMIS_MASS_IN,
     &                  PAQSO4,PH2SO4,PHNO3,PNH3,MH2SO4,MHNO3,MNH3,ITMATRIX
     
     	print*, ' SET', set
     	print*, ' NREP', nrep
     	print*, ' tstep', tstep
     	print*, ' tstop', tstop
    
     
          WRITE(*,'(A,2I8)')'Processing case N, SET = ', N, SET
        IF( WRITE_LOG ) WRITE(AUNIT1,'(//A,I5)')'NEW CASE: SET = ', SET
        IF( WRITE_LOG ) WRITE(AUNIT1,'(/A/A)')'SET,ICSET,ILAY,NREP,TSTEP,TSTOP,TK,RH,PRES,EMIS_MASS_IN,',
     &                                        'PAQSO4,PH2SO4,PHNO3,PNH3,MH2SO4,MHNO3,MNH3,ITMATRIX'
        IF( WRITE_LOG ) WRITE(AUNIT1,90005)    SET,ICSET,ILAY,NREP,TSTEP,TSTOP,TK,RH,PRES,EMIS_MASS_IN,
     &                                         PAQSO4,PH2SO4,PHNO3,PNH3,MH2SO4,MHNO3,MNH3,ITMATRIX
        IF( WRITE_LOG ) WRITE(34,'(/A,I5)')'SET = ', SET
        IF( WRITE_LOG ) WRITE(35,'(/A,I5)')'SET = ', SET
        IF( WRITE_LOG ) WRITE(41,'(/A,I5)')'SET = ', SET

        IF ( WRITE_OUT ) THEN
          WRITE(DATCONV,'(I2.2)') SET
          READ (DATCONV,'(A2  )') NSETCHAR
          IF( DISCRETE_EVAL_OPTION ) THEN
            FILENAME = 'set'//NSETCHAR//'disc.txt'
            OPEN(37,STATUS='REPLACE',FILE=TRIM(FILENAME))
            FILENAME = 'set'//NSETCHAR//'pdfi.txt'  
            OPEN(38,STATUS='REPLACE',FILE=TRIM(FILENAME))
            FILENAME = 'set'//NSETCHAR//'pdff.txt'
            OPEN(39,STATUS='REPLACE',FILE=TRIM(FILENAME))
          ENDIF
          FILENAME = 'set'//NSETCHAR//'pdfn.txt'
          OPEN(42,STATUS='REPLACE',FILE=TRIM(FILENAME))
          FILENAME = 'set'//NSETCHAR//'pdfm.txt'
          OPEN(43,STATUS='REPLACE',FILE=TRIM(FILENAME))
          DO I=1, NFILES
            FILENAME = 'set'//NSETCHAR//EXTCHAR(I)//'.txt'
            OPEN(OUTUNIT(I), STATUS='REPLACE',FILE=TRIM(FILENAME))
            IF ( EXTCHAR(I) .EQ. 'spcm' ) THEN
              !--------------------------------------------------------------------------------------------------------
              ! WRITE(UNIT2,*)'echo '//QUOTE//TRIM(FILENAME)//QUOTE
              ! WRITE(UNIT2,*)'diff ../group1/',TRIM(FILENAME),' ../group2/',TRIM(FILENAME)
              !--------------------------------------------------------------------------------------------------------
            ENDIF
          ENDDO
        ELSE
          IF( DISCRETE_EVAL_OPTION ) THEN
            OPEN(37,STATUS='REPLACE',FILE='disc.txt')
            OPEN(38,STATUS='REPLACE',FILE='pdfiniti.txt')
            OPEN(39,STATUS='REPLACE',FILE='pdffinal.txt')
          ENDIF
        ENDIF
 
        DT = TSTEP / 3600.0D+00                       ! convert [s] to [h]
        NIT =  TSTOP / DT

        EMIS_MASS(1:10) = EMIS_MASS_IN                ! mass emission rates         [ug/m^3/s]

        AQSO4RATE = PAQSO4                            ! incloud SO4 production rate [ug/m^3/s]
        ! GAS_IN( GAS_H2SO4 ) = MH2SO4 * 1.59517D-10  ! convert [molecules/cm^3] to [ugSO4/m^3]
        GAS_IN( GAS_H2SO4 ) = MH2SO4                  ! [ugSO4/m^3] directly input
        GAS_IN( GAS_HNO3  ) = MHNO3                   ! [ug/m^3]
        GAS_IN( GAS_NH3   ) = MNH3                    ! [ug/m^3]

        CALL INIT_AERO( ICSET, TK, PRES )           ! get initial aerosol concentrations
  
        IF ( WRITE_LOG ) THEN
          WRITE(AUNIT1,'(/A )') 'INITIAL AEROSOL CONCENTRATIONS'          
          WRITE(AUNIT1,'(/A/)') 'I, AERO_SPCS(I), AERO_IN(I) [#/m^3]'
          DO I=1, NAEROBOX 
            IF (AERO_SPCS(I)(1:4) .EQ. 'NUMB') WRITE(AUNIT1,90000) I, AERO_SPCS(I), AERO_IN(I)
          ENDDO
          WRITE(AUNIT1,'(/A/)') 'I, AERO_SPCS(I), AERO_IN(I) [ug/m^3]'
          DO I=1, NAEROBOX  
            IF (AERO_SPCS(I)(1:4) .EQ. 'MASS') WRITE(AUNIT1,90000) I, AERO_SPCS(I), AERO_IN(I)
          ENDDO
        ENDIF

        DO IREP=1, NREP

          AERO(:)  = AERO_IN(:)
          GAS(:)   = GAS_IN(:)   ! must come after call to SPCMASSES
          T = 0.0D+00
          DIAG_ACCUM(:,:) = 0.0D+00

          CALL SPCMASSES(AERO,GAS,SPCMASS)
          IF ( WRITE_OUT ) CALL OUT(AERO,GAS,NACT,SPCMASS,T,OUTUNIT)

          DO IT=1, NIT

            GAS( GAS_H2SO4 ) = GAS( GAS_H2SO4 ) + PH2SO4*TSTEP
            GAS( GAS_HNO3  ) = GAS( GAS_HNO3  ) + PHNO3 *TSTEP
            GAS( GAS_NH3   ) = GAS( GAS_NH3   ) + PNH3  *TSTEP

            IF ( CPU_STATS ) CALL CPUTIME( 1 )

            IF( DISCRETE_EVAL_OPTION ) CALL DISCRETE( ICSET, T, TSTEP )

            ! WRITE(*,*)'TOT_SULF (DRIVER) = ', SUM( AERO( SULF_MAP(:)) )
            ! WRITE(*,'(8F12.6)')AERO( SULF_MAP(:)) 

            IF( WRITE_LOG ) WRITE(40,'(F8.2)') T
            TSUBSTEP = TSTEP / REAL( ITMATRIX )
            DO I=1, ITMATRIX
              CALL MATRIX(AERO,GAS,EMIS_MASS,TSUBSTEP,TK,RH,PRES,AQSO4RATE,WUPDRAFT,DIAG,NACT) 
              DIAG_ACCUM(:,:) = DIAG_ACCUM(:,:) + DIAG(:,:) 
            ENDDO
            
            ! WRITE(*,*)'TOT_SULF (DRIVER) = ', SUM( AERO( SULF_MAP(:)) )
            ! WRITE(*,'(8F12.6)')AERO( SULF_MAP(:)) 

!           IF( WRITE_LOG ) WRITE(AUNIT1,'(/A,F12.6/)')'DIAM(1,1,10,1)*1.0D+06',DIAM(1,1,10,1)*1.0D+06

            T = T + DT
            ! IF ( T .GT. TEMIS_MAX .AND. SET .EQ. 15 ) EMIS_MASS(:) = 0.0D+00 

!----------------------------------------------------------------------------------------------------------------------
!           IF( DO_DIAG_WRITE ) THEN
!             WRITE(44,*)'TIME (h) = ', T
!             WRITE(44,*)'Diagnostics for number concentration'
!             DO I=1, NAEROBOX
!               IF( AERO_SPCS(I)(1:4) .NE. 'NUMB' ) CYCLE
!               WRITE(44,90007) I,AERO_SPCS(I),DIAG_ACCUM(1:7,I),SUM(DIAG_ACCUM(1:6,I))
!             ENDDO
!             WRITE(44,*)'Diagnostics for mass concentration'
!             DO I=1, NAEROBOX
!               IF( AERO_SPCS(I)(1:4) .EQ. 'NUMB' ) CYCLE
!               WRITE(44,90007) I,AERO_SPCS(I),DIAG_ACCUM(8:15,I),SUM(DIAG_ACCUM(8:14,I))
!             ENDDO
!           ENDIF
!----------------------------------------------------------------------------------------------------------------------

            CALL SPCMASSES(AERO,GAS,SPCMASS)
            IF ( WRITE_OUT ) CALL OUT(AERO,GAS,NACT,SPCMASS,T,OUTUNIT)

            IF( UPDATE_VDEP ) THEN
              DO I=1, NMODES
                WRITE(41,90006) T, I, MODE_NAME(I), 1.0D+06*DIAM(IXXX,IYYY,ILAY,I), VDDEP_AERO(IXXX,IYYY,I,:)
               !WRITE(41,90006) T, I, MODE_NAME(I),                        DGN0(I), VDDEP_AERO(IXXX,IYYY,I,:)
              ENDDO
              WRITE(41,*) '   '
            ENDIF

          ENDDO

          IF( DO_DIAG_WRITE ) THEN
            DO I=1, NAEROBOX
            DO J=1, NDIAG_AERO
              IF( ABS( DIAG_ACCUM(J,I) ) .LT. 1.0D-09 ) DIAG_ACCUM(J,I) = 0.0D+00 
            ENDDO
            ENDDO
            WRITE(44,*)'SET, TIME (h) = ', SET, T
            WRITE(44,*)'Diagnostics for number concentration'
            DO I=1, NAEROBOX
              IF( AERO_SPCS(I)(1:4) .NE. 'NUMB' ) CYCLE
              SUMTMP = 0.0D+00
              DO J=1,6
                SUMTMP = SUMTMP + DIAG_ACCUM(J,I)
              ENDDO
              WRITE(44,90007) I,AERO_SPCS(I),DIAG_ACCUM(1:7,I),SUMTMP
            ENDDO
            WRITE(44,*)'Diagnostics for mass concentration'
            DO I=1, NAEROBOX
              IF( AERO_SPCS(I)(1:4) .EQ. 'NUMB' ) CYCLE
              SUMTMP = 0.0D+00
              DO J=8,14
                SUMTMP = SUMTMP + DIAG_ACCUM(J,I)
              ENDDO
              WRITE(44,90007) I,AERO_SPCS(I),DIAG_ACCUM(8:15,I),SUMTMP
            ENDDO
          ENDIF

        ENDDO

        CALL PLOT_FINAL_PDFS(ICSET,42,43,AERO)
        IF( DISCRETE_EVAL_OPTION ) CALL DISCRETE_OUT(39,ICSET,T)

        DO I=1, NFILES
          CLOSE( OUTUNIT(I) )
        ENDDO
        CLOSE( 37 )
        CLOSE( 38 )
        CLOSE( 39 )
        CLOSE( 42 )
        CLOSE( 43 )

        IF ( CPU_STATS ) CALL CPUTIME( 1 )

        IF ( MINVAL( AERO ) .LT. 0.0D+00 .OR. MINVAL( GAS ) .LT. 0.0D+00 ) THEN
          WRITE(*,*)'Negative value in AERO array or GAS array: ', GAS(:)
          DO I=1, NAEROBOX
            WRITE(*,'(I5,D15.5,5X,A)')I,AERO(I),AERO_SPCS(I)
          ENDDO
          DO I=1, NGASES
            WRITE(*,'(I5,D15.5,5X,A)')I,GAS(I)
          ENDDO
          STOP
        ENDIF
        IF ( MAXVAL( AERO ) .GT. 1.0D+20 .OR. MAXVAL( GAS ) .GT. 1.0D+20 ) THEN
          WRITE(*,*)'Large value in AERO array or GAS array: ', GAS(:)
          DO I=1, NAEROBOX
            WRITE(*,'(I5,D15.5,5X,A)')I,AERO(I),AERO_SPCS(I)
          ENDDO
          DO I=1, NGASES
            WRITE(*,'(I5,D15.5,5X,A)')I,GAS(I)
          ENDDO
          STOP
        ENDIF

      ENDDO   ! DO N=1, NSETS
      
      CALL WRITE_DIAM_HISTOGRAM      

      IF ( CPU_STATS ) THEN
        CALL CPUTIME( -1 )
      ENDIF

      IF( DO_TEST_COAG  ) CALL TEST_COAG_COEF
      IF( DO_TEST_QUAD  ) CALL TEST_QUAD
      IF( DO_TEST_ACTV  ) CALL TEST_GETACTFRAC
      IF( DO_TEST_DEPV  ) CALL TEST_DEPV

      CLOSE(UNIT1)
      CLOSE(UNIT2)
      DO I=31, 41
        CLOSE(I)
      ENDDO
      CLOSE(44)
      CLOSE(AUNIT1)

      !-------------------------------------------------------------------------
      ! WRITE(*,*) NUMB_MAP(:)
      ! WRITE(*,*) SULF_MAP(:)
      ! WRITE(*,*) BCAR_MAP(:)
      ! WRITE(*,*) OCAR_MAP(:)
      ! WRITE(*,*) DUST_MAP(:)
      ! WRITE(*,*) SEAS_MAP(:)
      ! WRITE(*,*) MODE_NUMB_SEAS(:)
      ! WRITE(*,*) NM(:)
      ! DO I=1, NMODES
      !   WRITE(*,'(5A5)') NM_SPC_NAME(I,:)
      ! ENDDO
      !-------------------------------------------------------------------------

90000 FORMAT(I4,4X,A16,D15.5)
90001 FORMAT(I4,2D15.5,A20)
90002 FORMAT(8F14.7)
90003 FORMAT(14D13.5)
90004 FORMAT(14D13.5)
90005 FORMAT(4I3,4F8.2,F9.1,D12.4,/,7D12.4,I6)
90006 FORMAT(F10.4,I4,4X,A6,4X,3F13.8)
90007 FORMAT(I3,2X,A13,9D12.4)
      STOP
      END PROGRAM DRIVER


      SUBROUTINE OUT(AERO,GAS,NACT,SPCMASS,T,OUTUNIT)
      USE AERO_CONFIG
      USE AERO_SETUP  
      USE AERO_SUBS
      USE AMP_AEROSOL
      IMPLICIT NONE
      INTEGER :: I,J
      REAL(8) :: AERO(:)                   
      REAL(8) :: GAS(:)          
      REAL(8) :: NACT(:)          
      REAL(8) :: SPCMASS(:) ! total mass conc. of each model species [ug/m^3]
      REAL(8) :: T          ! model time [h]
      REAL(8) :: X(NMODES)            
      REAL(8) :: XN(NWEIGHTS),XM(NWEIGHTS)
      REAL(8) :: D(NWEIGHTS)            
      INTEGER :: OUTUNIT(:)
      REAL(8), PARAMETER :: MIN_NUMB          = 1.0D+02                       ! [ #/m^3]
      REAL(8), PARAMETER :: MASS_PER_PARTICLE = 5.236D-16                     ! [ug] for Dp=1nm and density=1g/cm^3
      REAL(8), PARAMETER :: MIN_MASS          = MIN_NUMB * MASS_PER_PARTICLE  ! [ug/m^3]
      INTEGER :: HALFN

      IF ( NMODES .EQ.  8 ) HALFN = 7   ! config 8
      IF ( NMODES .EQ. 11 ) HALFN = 7   ! config 7
      IF ( NMODES .EQ. 14 ) HALFN = 7   ! config 5, 6
      IF ( NMODES .EQ. 10 ) HALFN = 8   ! config 4
      IF ( NMODES .EQ. 13 ) HALFN = 8   ! config 3
      IF ( NMODES .EQ. 16 ) HALFN = 8   ! config 1, 2

      ! Total dry mass in each mode [ug/m^3].
	

      DO I=1, NMODES
        X(I) = 0.0D+00
        DO J=1, NAEROBOX
          IF ( AERO_SPCS(J)(1:8) .EQ. 'MASS_'//MODE_NAME(I) ) THEN
            !WRITE(*,*) I, AERO_SPCS(J),'MASS_'//MODE_NAME(I)
            WRITE(*,*) I,J,MODE_NAME(I),AERO_SPCS(J)
            X(I) = X(I) + AERO( J )
          ENDIF
        ENDDO
        X(I) = MAX( X(I), MIN_MASS )
      ENDDO
      WRITE(OUTUNIT( 1),90000) T,X(1:HALFN)
      WRITE(OUTUNIT( 2),90000) T,X(HALFN+1:NMODES)

      ! Number conc. for each mode [#/m^3].

      XN(:) = MAX( AERO( NUMB_MAP(:) ), MIN_NUMB ) * 1.0D-06  ! [#/m^3] to [#/cm^3]
      WRITE(OUTUNIT( 3),90000) T, XN(1:HALFN)
      WRITE(OUTUNIT( 4),90000) T, XN(HALFN+1:NWEIGHTS)

      ! Sulfate conc. for each mode [ug/m^3].

      XM(:) = MAX( AERO( SULF_MAP(:) ), MIN_MASS )
      WRITE(OUTUNIT( 5),90000) T, XM(1:HALFN)
      WRITE(OUTUNIT( 6),90000) T, XM(HALFN+1:NWEIGHTS),GAS( GAS_H2SO4 )

      ! BC, OC, DUST, SEAS conc. for each mode [ug/m^3].

      XM(:) = 0.0D+00
      XM(1:SIZE(BCAR_MAP)) = MAX( AERO( BCAR_MAP(:) ), MIN_MASS )
      WRITE(OUTUNIT( 7),90000) T, XM(1:SIZE(BCAR_MAP))

      XM(:) = 0.0D+00
      XM(1:SIZE(OCAR_MAP)) = MAX( AERO( OCAR_MAP(:) ), MIN_MASS )
      WRITE(OUTUNIT( 8),90000) T, XM(1:SIZE(OCAR_MAP))

      XM(:) = 0.0D+00
      XM(1:SIZE(DUST_MAP)) = MAX( AERO( DUST_MAP(:) ), MIN_MASS )
      WRITE(OUTUNIT( 9),90000) T, XM(1:SIZE(DUST_MAP))

      XM(:) = 0.0D+00
      XM(1:SIZE(SEAS_MAP)) = MAX( AERO( SEAS_MAP(:) ), MIN_MASS )
      WRITE(OUTUNIT(10),90000) T, XM(1:SIZE(SEAS_MAP))

      ! Diameter of average mass for each mode [um].

      D(:) = DP0(:)*1.0D+06
      DO I=1, NMODES
        ! WRITE(*,*) X(I),AERO( NUMB_MAP(I) )
        IF ( AERO( NUMB_MAP(I) ) .LE. 1.0D-05 ) CYCLE
        D(I) = 1.0D+06 * ( CONV_MASS_TO_DP * X(I)/AERO( NUMB_MAP(I) ) )**0.33333333333
      ENDDO
      WRITE(OUTUNIT(11),90000) T,D(1:HALFN)
      WRITE(OUTUNIT(12),90000) T,D(HALFN+1:NMODES)

      ! Total mass conc. of each chemical species: SO4,BC,OC,DUST,SEAS,NO3,NH4 [ug/m^3].

!     WRITE(OUTUNIT(13),90000) T, MAX( SPCMASS(:), MIN_MASS )
      WRITE(OUTUNIT(13),90000) T, MAX( SPCMASS(:), 0.0D+00  )
!     Gas concentrations
      WRITE(OUTUNIT(14),90000) T, MAX( GAS(:), 0.0D+00  )
!     WRITE(35,90000)          T, MAX( SPCMASS(:), MIN_MASS )
      WRITE(35,90000)          T, MAX( SPCMASS(:), 0.0D+00  )
      
      !CCN conc. for each mode [#/m^3].
      XM(:) = 0.0D+00
       print*, NACT
      XM(1:16) = MAX( NACT(:), MIN_NUMB) * 1.0D-06 ! [#/m^3] to [#/cm^3]
      XM(1:16) = MIN( XM(1:16), 1.d33) 
       print*, XM

!      XM(1:16) = MAX( XM(1:16) , MAX_NUMB ) 
!      WRITE(OUTUNIT( 15),90000) T, XM(1:16)
      WRITE(OUTUNIT( 15),90000) T,XM(1:HALFN)
      WRITE(OUTUNIT( 16),90000) T,XM(HALFN+1:NMODES)      
      
90000 FORMAT(F8.4,9D13.6)
      RETURN
      END SUBROUTINE OUT
      
      
      SUBROUTINE CPUTIME(ICPUFLAG)
      USE AERO_PARAM, ONLY: AUNIT1
!-------------------------------------------------------------------------------
!
!     Keeps CPU time statistics for selected sections of MATRIX.
!
!     CALL P2CPUTIM( 0) for initialization
!     CALL P2CPUTIM( I) for Section I
!     CALL P2CPUTIM(-1) for write of final statisitics to logical unit AUNIT1.
!
!-------------------------------------------------------------------------------
      IMPLICIT NONE
      INTEGER :: ICPUFLAG,I
      INTEGER, PARAMETER :: NCPUGRPS=15
      CHARACTER(LEN=10) :: CPUGRPNM(NCPUGRPS)
      CHARACTER(LEN=31) :: CPUSECTIONS(NCPUGRPS)
      REAL(8), SAVE :: TCPUGRP(NCPUGRPS),TCPUOLD,TCPU00,TCPU,TCPUTOT
      DATA CPUGRPNM /
     & 'SECTION  1', 'SECTION  2', 'SECTION  3',
     & 'SECTION  4', 'SECTION  5', 'SECTION  6',
     & 'SECTION  7', 'SECTION  8', 'SECTION  9',
     & 'SECTION 10', 'SECTION 11', 'SECTION 12',
     & 'SECTION 13', 'SECTION 14', 'SECTION 15'/
      DATA CPUSECTIONS /
     & 'DRIVER                         ',
     & 'UP TO CALL AERO_THERMO         ',
     & 'AERO_THERMO + HYSTERESIS       ',
     & 'PARTICLE MASSES                ',
     & 'UPDATE COAGULATION COEFFICIENTS',
     & 'DN/DT:  EMIS + NPF TERMS       ',
     & 'DN/DT:  COAG TERMS             ',
     & 'DM/DT:  EMIS TERMS             ',
     & 'DM/DT:  COAG TERMS             ',
     & 'DM/DT:  NUCL+CLOUD+GROWTH TERMS',
     & 'DN/DT:  SOLVER                 ',
     & 'DM/DT:  SOLVER                 ',
     & 'AERO_THERMO + HYSTERESIS       ',
     & 'INSOLUBLE-SOLUBLE MODE TRANSFER',
     & '(EMPTY)                        '/
 
      CALL CPUSECND(TCPU)

      IF (ICPUFLAG.EQ.0) THEN  ! First call - initialize.

        TCPU00 = TCPU
        TCPUGRP(:) = 0.0D+00

      ELSEIF(ICPUFLAG.GT.0) THEN  ! If flag is positive, store execution time for
                                  ! the cpu timing section designated by it.
        I = MIN(ICPUFLAG,NCPUGRPS)
        TCPUGRP(I) = TCPUGRP(I) + (TCPU - TCPUOLD)

      ELSE  ! If flag is negative, write total execution time for entire program
            ! for each of the cpu timing sections defined by CPUGRPNM.

        TCPU = TCPU - TCPUGRP(1)    ! remove driver CPU time before stats.
        TCPUGRP(1) = 0.0d+00        ! remove driver CPU time before stats.
        TCPUTOT = MAX( (TCPU - TCPU00), 1.0D-10 )
        WRITE(AUNIT1,9100) TCPUTOT, (TCPUTOT/60.0D+00)
        DO I = 1, NCPUGRPS
          WRITE(AUNIT1,9110) CPUGRPNM(I),TCPUGRP(I),100.0D+00*(TCPUGRP(I)/TCPUTOT),
     &                       CPUSECTIONS(I)
        ENDDO
        WRITE(AUNIT1,'(A)')'  '

      ENDIF

      TCPUOLD = TCPU

 9100 FORMAT(/'CPU TIME STATISTICS:' //
     &     '    TOTAL TIME (SECS, MINS)', 2(1X,F12.4) //
     &     '    SECTION           SECS      % OF TOTAL    ',
     &     '    CODE SECTION DESCRIPTION' / )
 9110 FORMAT(4X,A,2(1X,F12.4),10X,A31)

      RETURN
      END SUBROUTINE CPUTIME


      SUBROUTINE CPUSECND(TCPU)
!-------------------------------------------------------------------------------
!     Returns elapsed program execution time (s) since beginning of program
!     execution, in variable TCPU.
!-------------------------------------------------------------------------------
      IMPLICIT NONE
      REAL(8) :: TCPU
!-------------------------------------------------------------------------------
!     Uncomment this line if the function ETIME is not available.
!-------------------------------------------------------------------------------
      TCPU = 0.0D+00
!-------------------------------------------------------------------------------
!     Uncomment these three lines if the function ETIME is available.
!-------------------------------------------------------------------------------
!     REAL    :: TCPUSUN(2)
!     CALL ETIME(TCPUSUN)
!     TCPU = TCPUSUN(1) + TCPUSUN(2)
!-------------------------------------------------------------------------------
      RETURN
      END SUBROUTINE CPUSECND

c -----------------------------------------------------------------
      SUBROUTINE SPCMASSES(AERO,GAS,SPCMASS)
!----------------------------------------------------------------------------------------------------------------------
!     Routine to calculate the total mass concentration of each model species:
!     SULF, BCAR, OCAR, DUST, SEAS, NO3, NH4. Aerosol water is not treated. 
!----------------------------------------------------------------------------------------------------------------------
      USE AERO_SETUP, ONLY: SULF_MAP, BCAR_MAP, OCAR_MAP, DUST_MAP, SEAS_MAP
      USE AERO_PARAM
      USE AERO_CONFIG
      IMPLICIT NONE
      REAL(8) :: AERO(NAEROBOX)
      REAL(8) :: GAS(NGASES)    
      REAL(8) :: SPCMASS(NMASS_SPCS+2)
      SPCMASS(1) = SUM( AERO( SULF_MAP(:) ) ) + GAS( GAS_H2SO4 )
      SPCMASS(2) = SUM( AERO( BCAR_MAP(:) ) )
      SPCMASS(3) = SUM( AERO( OCAR_MAP(:) ) )
      SPCMASS(4) = SUM( AERO( DUST_MAP(:) ) )
      SPCMASS(5) = SUM( AERO( SEAS_MAP(:) ) )
      SPCMASS(6) = AERO( MASS_NO3 )           + GAS( GAS_HNO3 )
      SPCMASS(7) = AERO( MASS_NH4 )           + GAS( GAS_NH3  )

 
      RETURN
      END SUBROUTINE SPCMASSES
